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INTRODUCTION 


This  final  report  is  for  work  carried  out  under  Grant  No.  F49620-93-T0143 
during  the  18-month  period  from  1  April,  1993  to  30  September,  1994. 

The  principal  objective  of  the  research  was  to  make  use  of  a  parallel  com¬ 
puter  code  recently  developed  by  our  group,  for  Bird’s  direct  simulation  Monte 
Carlo  (DSMC)  method,  and  to  gain  experience  in  simulating  various  rarefied  three- 
dimensional  flows  about  bodies  having  reasonably  complex  geometries,  for  example, 
bodies  based  on  geometries  associated  with  current  practical  applications. 

The  computational  facilities  needed  to  carry  out  the  work  were  made  available 
as  a  result  of  our  close  collaboration  with  several  researchers  at  NASA-Ames  Research 
Center,  in  particular  W.  J.  Feiereisen  (Computational  Aerosciences  Project  Manager) 
and  T.  A.  Edwards  (Chief,  Aerothermodyncimics  Branch).  In  the  18-month  period 
of  the  grant  this  gave  our  students  access  to  the  use  of  supercomputers  such  as  the 
Intel  iPSC/860  Gamma,  Intel  iPSC/860  Delta,  the  Intel  Paragon  XP/S-15  and  the 
Cray  Y-MP  C90.  In  early  stages  of  their  work,  the  students  would  either  make 
use  of  an  HP  9000/730  workstation  in  our  laboratory  or  access  the  supercomputers 
from  microcomputers  in  our  laboratory.  In  later  stages,  they  had  the  opportunity  to 
conduct  their  research  on  workstations  at  NASA-Ames. 

The  work  carried  out  in  this  period  was  closely  tied  to  earlier  contributions 
of  several  students  who  participated  in  developing  the  version  of  the  code  used  in 
these  studies.  Before  the  start  of  this  investigation,  Jeffrey  D.  McDonald  had  finished 
a  multiple  species  version  of  a  particle  code  which  he  developed  to  take  advantage 
of  the  vector  capabilities  of  the  Cray  family  of  supercomputers  (Cray-2,  Cray  Y- 
MP),  as  well  as  to  take  advantage  of  certain  algorithmic  improvements  introduced 
by  our  group  [1,2].  The  most  recent  version  of  this  serial  code  is  called  PSim2. 
McDonald  then  worked  on  a  parallel  version  of  the  same  code  for  use  on  the  Intel 
family  of  parallel  supercomputers.  The  most  recent  version  of  this  code,  with  very 
important  contributions  from  Michael  A.  Fallavollita  [3],  is  called  PSim4.  Some 
additional  development  of  the  parallel  code  was  then  continued  by  Douglas  C.  Dahlby, 
to  introduce  changes  that  would  make  it  easier  to  use  in  the  study  of  problems  having 
important  practical  aspects. 
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PARALLEL  CODE  APPLICATIONS 


In  an  effort  to  gain  experience  in  the  application  of  our  parallel  computer  code  to 
problems  that  contain  practical  elements  with  moderately  complex  bodies,  we  turned 
our  attention  to  the  study  of  a  body  shape  that  corresponds  to  a  wing-body,  single- 
stage  to  orbit  (SSTO)  vehicle,  where  we  have  attempted  to  reproduce  the  general 
shape  of  a  proposed  design  but  have  used  our  own  values  for  the  different  parameters 
defining  the  vehicle  geometry.  Figure  1  shows  the  generic  vehicle  shape  we  have  been 
considering.  In  this  case  the  Mach  number  was  set  at  10,  the  angle  of  attack  at  20 
degrees,  and  the  Knudsen  number  at  0.05,  based  on  the  length  of  the  body.  The 
display  shows  the  density  distribution  in  the  central  plane  of  the  body  and  in  a  single 
cross  section  of  the  flow  just  downstream  of  the  body.  One  of  the  many  quantities 
that  can  be  obtained  from  the  simulation  is  the  heat-flux  to  the  body  surface,  which 
is  shown  in  Fig.  2  for  the  vehicle  at  zero  angle  of  attack.  The  flow  field  and  body 
were  resolved  by  use  of  cubic  cells,  with  120  in  both  the  stream-wise  and  vertical 
directions,  and  with  75  cells  along  the  semispan  direction,  for  a  total  of  1.1  million 
cells.  Because  we  did  not  feel  it  was  appropriate  to  invest  a  large  fraction  of  our  time 
in  the  detailed  design  of  a  body,  certain  details  such  as  the  inclusion  of  wing-body 
fillets  were  ignored.  After  the  file  was  created  that  defined  the  body,  the  generation 
of  the  three-dimensional  displays  shown  in  the  two  figures  was  made  possible  by  use 
of  a  graphics  software  package  called  MATLAB  (Math  Works  Inc.)  which  was  run  on 
our  HP  9000/730  UNIX  workstation;  color  hard-copy  prints  were  made  using  an  HP 
1200C  postscript  color  ink-jet  printer.  In  both  cases,  the  acquisition  and  maintenance 
costs  proved  to  be  nominal  and  ideal  for  a  university  environment. 

Since  the  SSTO  body  was  the  first  wing-body  we  studied,  we  were  interested 
in  determining  whether  a  vortical  flow  would  be  seen  in  the  DSMC  method  owing 
to  lift  on  the  wings.  Figures  3  and  4  show  the  velocity  fields  generated  at  cross 
sections  located  at  approximately  60%  and  90%  of  the  wing  root  chord,  respectively. 
The  simulation  was  carried  out  for  the  Ccise  of  a  subsonic  vehicle  Mach  number, 
where  vortical  structures  would  be  expected  to  be  clearly  present.  The  case  shown 
corresponds  to  an  angle  of  attack  of  20  degrees  and  a  Mach  number  of  0.8;  all  other 
parameters  were  the  same  as  those  employed  for  Fig.  1.  As  can  be  seen,  vortex 
structures  are  not  only  observed  for  the  entire  flow  but  one  can  also  see  in  the  second 
downstream  cross-section  that  certain  portions  of  the  vortex  flow  can  be  associated 
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with  the  main  wing  section  and  certain  portions  with  the  end  winglet.  Because  of 
the  necessity  of  using  finite  lengths  in  plotting  the  velocity  vectors,  visual  correction 
is  often  needed  for  vectors  located  near  solid  boundaries,  if  the  proper  boundary 
conditions  are  to  be  seen.  In  view  of  the  rarefied  conditions  of  the  run,  the  results  are 
very  fascinating  and  add  greatly  to  ones  intuitive  understanding  of  fluid  mechanics. 

For  problems  involving  steady  flow,  time  averaging  is  used  in  the  DSMC  method 
to  improve  the  simulated  results,  but  time  averaging  can  only  be  initiated  after  the 
steady-state  condition  has  been  reached.  For  some  problems,  the  transient  period  can 
be  a  significant  fraction  of  the  corresponding  averaging  time.  If  one  first  wishes  to 
explore  a  particular  problem  without  investing  a  large  amount  of  time  in  averaging, 
and  later  wishes  to  refine  the  results  by  increasing  the  averaging  time,  or  to  use 
the  results  as  a  starting  state  for  a  similar  problem,  then  the  ability  to  restart  the 
simulation  is  clearly  needed.  In  the  DSMC  method,  the  amount  of  data  that  must  be 
saved  and  then  loaded  for  a  restart  on  a  large  parallel  machine  is  truly  huge  and  quite 
impractical  to  consider.  However,  there  is  reason  to  believe  that  the  normal  output 
from  a  typical  run  can  be  used  to  re-create  an  initial  run  state  that  is  close  enough 
to  the  final  steady  state  that  the  corresponding  transient  period  would  be  negligibly 
small.  This  approach  has  been  implemented  and  found  to  be  quite  practical  and  very 
effective. 

Figure  5  shows  the  result  of  a  test  of  the  restart  capability  for  a  problem  where 
the  Knudsen  layer  is  the  critical  element  in  the  flow.  When  the  mean  free  path 
length  is  relatively  large  in  a  Prandtl-Meyer  flow,  a  Knudsen  layer  develops  on  the 
down-stream  wall  bounding  the  expanding  two-dimensional  supersonic  flow.  The 
Knudsen  layer  is  fully  developed  only  for  the  steady  state.  One  highly  resolved  run  was 
conducted  to  serve  as  the  reference  solution  and  then  the  rms  temperature  difference 
was  found  for  two  other  runs;  one  starting  from  a  conventional  initial  free-streaming 
state  and  the  other  starting  from  a  state  created  from  the  reference  solution,  where  a 
Maxwellian  velocity  distribution  function  is  initially  employed  for  all  cells.  The  figure 
shows  that  the  reconstructed  initial  state  approaches  the  final  steady  state  much  more 
quickly  than  the  start-up  used  in  the  conventional  method.  A  second  restart  test  was 
conducted  for  the  case  where  a  shock  wave  became  the  critical  element  in  the  flow. 
In  this  test,  a  flow  past  a  circular  cylinder  was  studied.  Figure  6  shows  the  results 
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obtained  and  it  is  clear  that  in  this  case  a  considerable  number  of  time  steps  are  saved 
by  the  restart  procedure  in  reducing  the  transient  period. 

Because  the  restart  procedure  used  constructs  an  Euler  flow  field,  the  greatest 
error  is  clearly  introduced  in  shock  layers.  This  can  be  seen  in  Fig.  7,  where  the 
average  rms  density  error  between  the  restart  solution  and  the  reference  solution  is 
displayed  for  the  period  50  to  100  time  steps  after  restart.  The  greatest  difference  is 
clearly  seen  in  the  shock  wave  as  it  attempts  to  adjust  itself  to  the  correct  form,  while 
the  corresponding  error  in  the  boundary  layer  on  the  cylinder  is  not  seen  as  clearly 
on  the  scale  displayed  in  the  figure. 

The  smallest  volume  of  physical  space  used  in  the  simulation  is  a  cell,  usually 
cubic  in  shape.  Of  course  these  should  be  made  as  small  as  possible  for  good  resolution 
of  both  the  body  and  the  flow  field.  However,  in  certain  regions  where  the  density 
is  low,  there  may  be  far  too  few  particles  per  cell  to  correctly  model  the  collision 
mechanics,  and  thus  get  accurate  results,  even  when  using  lengthy  time  averaging.  In 
those  regions  it  would  be  useful  to  collect  several  adjacent  cells  into  a  single  macrocell 
while  carrying  out  the  simulation.  A  test  study  was  conducted  where  the  flow  ahead 
of  and  within  the  bow  shock  wave  generated  by  a  circular  cylinder  was  intentionally 
set  at  a  level  where  too  few  particles  per  cell  were  present  to  produce  accurate  results 
in  those  regions  when  using  single  cells.  A  second  test  study  was  then  carried  out  at 
exactly  the  same  conditions  except  macrocells  were  introduced  in  those  regions  where 
the  particle  count  was  too  low.  Both  of  these  runs  were  compared  with  a  reference 
run  carried  out  at  a  much  higher  particle  number  density. 

The  first  test  study  conducted  used  a  particle  density  of  1  particle  per  cell  in 
the  upstream  flow;  and  the  reference  run  was  carried  out  at  a  density  of  16  particle 
per  cell  in  the  upstream  flow.  These  results  are  shown  in  Fig.  8a  where  the  reference 
solution  is  given  by  the  solid  curves  and  the  test  study  is  given  by  the  dashed  curves. 
One  can  see  that  having  too  few  particles  per  cell  shifts  the  position  of  the  shock  wave 
toward  the  high  density  side  and  causes  the  shock  wave  to  become  too  thin,  leading 
to  error  in  the  prediction  of  the  flow  field.  The  second  test  study  was  carried  out  with 
the  same  particle  density  of  1  particle  per  cell  in  the  upstream  flow,  but  macrocells 
were  introduced  in  the  low-density  regions  to  correct  the  collision  mechanics  in  the 
simulation.  These  results  are  shown  by  the  dashed  curves  in  Fig.  8b  and  they  are  seen 
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to  compare  very  well  with  the  reference  solution,  establishing  the  fact  that  macrocells 
can  be  used  in  an  effective  way  to  improve  the  predictions  of  a  DSMC  simulation. 

In  this  comparison,  the  total  computation  time  was  held  fixed  by  increasing  the 
averaging  time  for  the  runs  having  the  low  particle  number  density.  The  important 
conclusion  drawn  is  that  when  computer  memory  limits  one  to  the  use  of  1  particle 
per  cell  in  a  given  simulation  then  the  introduction  of  macrocells  allows  one  to  obtain 
results  that  are  equivalent  in  quality  to  a  run  employing  16  times  as  many  particles 
per  cell.  That  is,  a  factor  of  16  reduction  in  the  size  of  computer  memory  required 
clearly  has  a  large  impact  on  the  selection  of  problems  that  can  be  run  for  a  given 
computer. 

Before  one  can  confidently  predict  the  optimum  operating  conditions  for  a  given 
simulation,  one  must  thoroughly  explore  parameter  space  for  a  wide  variety  of  con¬ 
ditions.  One  question  that  can  be  raised  is  how  important  is  the  trade-off  between 
the  improved  spatial  resolution  gotten  by  reducing  the  size  of  each  cell  versus  the  ef¬ 
fect  of  the  resultant  reduction  in  the  number  of  particles  per  cell,  assuming  the  total 
number  of  particles  and  thus  the  computational  effort  is  held  fixed?  Figure  9  gives 
a  summary  of  such  a  study,  where  the  problem  represents  the  simulation  of  three- 
dimensional  supersonic  flow  past  a  sphere,  and  the  quantity  studied  was  the  rms 
temperature  difference  with  respect  to  a  reference  solution.  The  size  of  the  sphere 
relative  to  the  size  of  the  wind  tunnel  was  held  fixed,  while  the  number  of  cells  making 
up  the  wind  tunnel  was  increased  by  factors  of  two,  from  a  resolution  of  8^  cells  to  64® 
cells.  When  a  total  of  1.05  million  particles  are  used.  Fig.  9a  showns  that  the  solution 
continually  improves  as  the  cell  size  is  decrecised,  indicating  that  improved  spatial 
resolution  is  the  controlling  factor,  even  though  an  average  as  low  as  approximately 
4  particles  per  cell  is  reached.  For  a  lesser  total  amount  of  particles,  the  trade-off  can 
be  clearly  seen,  with  an  optimum  occurring  at  32®  cells.  When  the  same  data  are 
displayed  as  in  Fig.  9b,  one  is  able  to  study  the  effect  of  increasing  the  total  number 
of  particles,  assuming  a  fixed  spatial  resolution  or  cell  size.  This  study  was  conducted 
to  explore  a  concept,  while  working  within  the  bounds  of  simulation  capability  and 
acceptable  run-time  cost  for  a  particular  problem,  rather  than  to  produce  results  that 
can  be  used  to  draw  general  conclusions.  A  full  report  on  the  above  material  together 
with  further  refinements  and  applications  of  our  parallel  code  will  appear  in  a  Ph.D. 
thesis  by  Douglas  Dahlby  within  the  next  year. 
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HYBRID  MODEL 


The  DSMC  method  becomes  computationally  inefficient  when  the  conditions  of 
the  flow  approach  the  near  continuum,  i.e.,  when  the  Navier-Stokes  equations  become 
valid.  In  this  limit  the  required  cell  size  becomes  small  and  the  number  of  particles 
required  by  the  simulation  becomes  large.  For  flow  conditions  where  the  stagnation 
enthalpy  is  very  high  and  the  body  temperature  is  low,  the  boundary  layer  may 
become  a  high-density,  near-continuum  region  that  is  best  handled  by  the  Navier- 
Stokes  equations.  Therefore,  it  is  most  desirable  to  develop  a  procedure  whereby  the 
DSMC  method  and  a  continuum  method  can  be  employed  in  their  respective  regions 
of  validity  and  then  joined  in  a  suitable  fashion  at  some  interface,  i.e.,  the  creation 
of  a  hybrid  scheme.  This  need  has  been  previously  considered  for  certain  special 
problems,  especially  in  the  study  of  plume  flows  in  space  where  the  gas  density  in  a 
nozzle  is  high  and  the  density  in  the  outer  plume  is  very  low,  clearly  suggesting  that 
the  two  approaches  be  joined  by  some  suitable  scheme. 

The  central  issue  in  introducing  a  hybrid  scheme  is  to  identify  the  proper  con¬ 
ditions  to  be  used  at  the  interface.  For  the  nozzle  problem,  where  the  flow  in  the 
nozzle  is  nearly  an  Euler  flow,  the  matching  conditions  are  not  critical  and  probably 
can  be  handled  by  assuming  a  Maxwellian  velocity  distribution  function,  to  model 
the  one-sided  fluxes  emanating  from  the  Euler  flow.  However,  in  flows  where  viscous 
and  heat  conduction  effects  dominate,  such  as  in  a  boundary  layer,  the  matching 
conditions  are  far  more  critical,  not  presently  known,  and  have  not  been  adaquately 
studied.  One  of  the  principal  problems  that  arises  in  this  case  is  that  the  physics 
associated  with  the  DSMC  method  leads  to  a  clear  definition  for  the  one-sided  fluxes 
(mass,  momentum,  and  energy)  passing  through  a  given  interface.  On  the  other  hand, 
the  Navier-Stokes  equations  provide  a  clear  definition  for  the  total  fluxes  (the  sum 
of  the  left-  and  right- directional  values).  Although  the  differences  between  the  two 
properly  define  the  fluxes  entering  the  DSMC  zone,  one  does  not  have  access  to  the 
proper  velocity  distribution  function  for  these  returning  particles,  and  this  is  needed 
if  the  DSMC  zone  is  to  blend  smoothly  with  the  Navier-Stokes  zone.  Additionally, 
if  one  were  to  make  use  of  one  of  the  flux  vector  splitting  schemes  in  computational 
fluid  dynamics  (CFD),  such  as  the  Steger- Warming  scheme,  the  split  fluxes  do  not 
have  the  same  physical  meaning  as  the  one-sided  fluxes  in  the  DSMC  method,  and 
therefore,  lead  to  new  complications. 
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Our  objective  in  studying  the  problem  was  to  identify  the  proper  general  con¬ 
ditions  to  be  used  at  the  interface  between  a  DSMC  simulation  and  a  suitable  CFD 
calculation.  The  initial  approach  taken  was  to  study  several  simple  problems  where 
viscous  and  heat  conduction  effects  are  clearly  important  and  then  explore  different 
matching  conditions  at  the  interface,  starting  with  the  simplest  concepts  and  then 
evolving  towards  the  use  of  more  complex  concepts,  as  required. 

The  first  model  problem  considered  was  a  strict  one-dimensional  nonsteady  flow 
in  a  finite  tube  with  closed  ends.  One  half  of  the  tube  was  initially  filled  with  an 
equilibrium  gas  and  the  other  half  remained  empty.  This  corresponds  to  the  shock- 
tube  problem  when  no  test  gas  is  present,  only  a  driver  gas.  When  the  driver  gas 
is  suddenly  released  by  the  bursting  diaphragm,  an  expansion  fan  forms  and  the 
accelerated  gas  rushes  towards  its  corresponding  closed  end.  On  striking  the  close 
end,  a  reflected  shock  wave  forms  and  propagates  back  into  the  onrushing  flow  created 
by  the  expansion  process.  The  hybrid  concept  is  studied  by  using  a  CFD  scheme  to 
handle  the  zone  represented  by  the  initial  driver  gas  and  the  DSMC  method  to  handle 
the  zone  represented  by  the  initial  vacuum.  A  comparison  is  made  by  employing  a 
reference  solution  obtained  by  using  the  DSMC  method  to  simulate  the  entire  flow. 

Results  of  such  a  study  are  presented  in  Fig.  10,  where  the  vertical  dash-dot 
line  identifies  the  location  of  the  fixed  interface  station  used,  with  the  CFD  zone  on 
the  left  side  and  the  DSMC  zone  on  the  right  side  of  each  plot.  Because  a  CFD 
solution  is  not  able  to  handle  regions  of  zero  density,  the  interface  station  had  to 
be  placed  slightly  to  the  left  of  the  diaphragm  station  (located  at  the  midpoint  of 
the  tube).  The  dashed  curve  gives  the  reference  DSMC  solution  and  the  solid  curve 
gives  the  hybrid  solution,  with  the  CFD  result  to  the  left  of  the  vertical  dash-dot 
line  and  the  DSMC  result  to  the  right  of  the  vertical  line.  The  two  simulations  were 
advanced  until  the  reflected  shock  wave  just  reached  the  interface  station  (580  time 
steps).  This  is  the  time  for  which  the  comparisons  shown  in  the  figure  are  made. 
For  these  initial  studies,  for  which  we  set  Kn  =  0.5  (cell  Knudsen  number),  CFL 
=  0.2  and  tube  length  =  100  cells,  we  assumed  that  the  conditions  at  the  interface 
could  be  approximated  by  using  a  Maxwellian  velocity  distribution  function.  Until 
the  reflected  shock  wave  arrives  at  the  interface  station,  the  assumption  proves  to  be 
a  good  one  since  the  expansion  fan  is  essentially  an  isentropic  process.  This  is  borne 
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out  by  the  data  in  the  figure  as  the  hybrid  solution  matches  the  reference  solution 
quite  well  for  all  variables  displayed. 

The  same  comparison  is  shown  in  Fig.  11  for  a  slightly  later  time  (620  time 
steps),  a  time  when  the  high  density  side  of  the  reflected  shock  wave  reaches  the 
interface  station.  At  this  point  the  simple  matching  conditions  used  are  clearly  in 
error  because  of  the  differences  seen  in  the  temperature  and  density  variables  at 
the  interface  station.  Similar  differences  are  seen  at  later  times,  such  as  when  the 
shock  wave  reaches  the  closed  end  at  the  left-hand  side,  as  well  as  when  it  passes  the 
interface  station  the  second  time.  This  test  proves  to  be  quite  useful  because  of  the 
otherwise  good  agreement  found  between  the  two  simulation  schemes;  it  is  only  at 
specific  locations  like  the  interface  station  that  the  error  is  highlighted. 

The  impulsively  started  flat  plate  represents  another  problem  for  which  useful 
testing  can  be  carried  out.  To  properly  handle  the  Knudsen  layer  that  always  develops 
close  to  a  material  surface,  we  must  place  the  DSMC  zone  next  to  the  plate  surface, 
while  the  CFD  zone  is  set  up  to  handle  the  remaining  outer  flow.  Selection  of  the 
interface  position  may  seem  to  be  arbitrary,  but  in  actual  fact  it  is  not  entirely 
arbitrary.  A  CFD  solution  depends  on  the  Navier-Stokes  equations  and  they  are 
valid  only  when  the  value  of  the  ratio  of  the  shearing  stress  to  the  local  pressure  r/p 
is  small.  Because  this  quantity  starts  as  a  very  large  number  next  to  the  plate  surface 
and  then  decays  as  it  diffuses  out  from  the  plate,  one  must  place  the  interface  station 
at  a  location  where  this  ratio  is  always  small.  If  not,  the  conditions  used  for  interfacing 
become  critical.  Such  a  test  is  shown  in  Fig.  12,  where  the  interface  station  is  located 
at  a  vertical  distance  of  0.2  mean  free  path  lengths  above  the  plate.  The  same  type  of 
comparison  with  a  hybrid  solution  is  made  as  that  considered  in  developing  Figs.  10 
and  11,  and  the  same  notation  for  each  curve  is  used,  except  in  this  case  we  are 
also  able  to  display  the  known  solution  for  the  incompressible  viscous  fluid,  shown 
as  a  dotted  curve.  The  observation  made  in  this  case  is  that  as  long  as  the  ratio 
r/p  is  smaller  than  0.1  beyond  the  interface  location,  the  hybrid  solution  compares 
favorably  with  the  reference  DSMC  result.  However,  when  the  ratio  becomes  equal 
to  0.1  at  the  interface  position,  the  hybrid  solution  and  the  reference  solution  show  a 
very  significant  difference,  as  seen  in  the  velocity  profile. 
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The  two  cases  studied  clearly  show  that  for  a  flow  that  is  dominated  by  vis¬ 
cous  stress  and  heat  conduction  the  interfacing  problem  requires  the  development  of 
an  approach  that  goes  beyond  the  simple  use  of  a  Maxwellian  velocity  distribution 
function  to  model  the  transfer  of  particles  from  a  CFD  solution  to  a  DSMC  solution. 
M^hat  is  needed  is  a  numerically  simple  scheme  that  can  account  for  the  translational 
nonequilibrium  of  the  gas  by  predicting  the  proper  number  of  particles  for  each  veloc¬ 
ity  interval,  as  particles  are  transferred  to  a  DSMC  solution.  This  work  was  carried 
out  by  Tawei  Lou  and  a  complete  report  will  be  given  in  his  Ph.D.  thesis,  which 
should  appear  in  the  next  year. 


COUETTE  FLOW:  A  REFERENCE  SOLUTION 

In  an  effort  to  better  understand  the  limitations  on  the  DSMC  method  when 
the  mean  free  path  length  becomes  small,  and  to  gain  the  necessary  information  for 
interfacing  a  DSMC  solution  with  a  continuum  solution,  we  turned  our  attention  to 
the  study  of  Couette  flow.  This  work  is  fully  described  in  a  Ph.D.  thesis  by  T.  Denery 
[4],  so  only  a  very  brief  description  of  the  results  obtained  will  be  given  here.  The 
interest  in  this  particular  problem  was  based  on  several  factors:  1)  Following  the 
recent  work  of  Baganoff  [5],  it  has  become  possible  to  extend  the  original  analytical 
work  of  Lees  and  Liu,  which  was  limited  to  the  case  of  Maxwell  molecules,  to  the  more 
general  case  of  inverse-power  molecules,  including  the  hard-sphere  limit;  2)  because 
of  our  previous  work  in  developing  a  DSMC  code  for  use  on  a  parallel  computer  [3], 
we  were  in  a  good  position  to  carry  out  a  large  number  of  simulations  of  the  flow 
for  comparison  purposes;  and  3)  The  DSMC  method  is  most  securely  founded  for  the 
case  of  hard-sphere  molecules,  a  case  for  which  a  theoretical  solution  became  available 
through  our  recent  work  [5],  and  therefore  it  has  become  possible  for  the  first  time  to 
carry  out  a  careful  and  proper  comparison  between  DSMC  and  theory. 

The  analytical  study  has  also  shed  light  on  a  very  important  question  in  the 
application  of  the  moment  method  in  kinetic  theory.  It  is  well  known  that  in  forming 
moments  of  the  Boltzmann  equation,  when  the  order  of  the  moment  is  higher  than 
the  orders  of  the  collisional  invariants,  that  the  corresponding  integration  required 
to  treat  the  collision  integral  cannot  be  carried  out  in  closed  form,  unless  one  as¬ 
sumes  an  interaction  potential  known  as  the  Maxwell  molecule.  A  common  practice 
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employed  when  considering  a  general  inverse-power  interaction  potential  has  been  to 
assume  that  one  can  used  the  closed  form  results  found  for  the  Maxwell  molecule 
and  simply  replace  the  linear  temperature  dependence  of  the  transport  coefficients 
for  the  Maxwell  molecule  by  the  particular  temperature  dependence  associated  with 
the  molecular  law  of  interest,  for  example,  square  root  for  the  case  of  the  hard-sphere 
molecule.  No  sound  theoretical  justification  exists  for  this  replacement,  but  it  has 
been  frequently  used  as  a  reasonable  approximation,  especially  in  the  analysis  of 
shock- wave  structure  and  its  comparison  with  experiment. 

Once  we  were  able  to  developed  an  analytic  solution  for  Couette  flow,  for  a 
general  inverse-power  molecule,  we  were  in  a  strong  position  to  test  the  above  pro¬ 
cedure.  Figure  13  shows  a  comparison  for  a  plate  Mach  number  of  0.5  and  for  two 
ratios  of  the  plate  temperature,  one  near  unity  and  the  other  fairly  extreme.  Both 
stress  and  heat  flux  are  shown,  normalized  by  their  respective  values  for  free-molecule 
flow,  with  the  Navier-Stokes  result  in  each  case  given  by  the  dotted  straight  line.  The 
comparison  of  interest  is  between  the  dashed  curve  and  the  solid  curve  in  each  case, 
which  are  seen  to  nearly  overlap;  and  because  the  hard-sphere  model  is  furthest  re¬ 
moved  from  the  Maxwell  molecule,  it  was  chosen  for  the  comparison.  The  solid  curve 
represents  the  theoretical  prediction  of  our  analytical  solution  for  hard  spheres,  while 
the  dashed  curve  gives  the  prediction  bcised  on  the  replacement  procedure  discussed 
above.  In  Fig.  14  the  same  comparison  is  made  except  the  plate  Mach  number  was 
raised  to  2.0.  Because  of  the  nearly  perfect  agreement  and  because  of  the  fact  that 
the  computational  effort  required  by  the  full  theory  is  over  10  times  greater  than  that 
required  for  the  replacement  procedure,  it  is  clear  that  the  common  practice  is  fully 
justified  for  Couette  flow;  and  it  lends  strong  support  to  the  belief  that  it  is  also  a 
very  reasonable  approximation  for  more  general  flows. 

One  of  the  most  remarkable  results  found  from  our  theoretical  work  for  Couette 
flow  is  exhibited  in  Fig.  15.  In  the  Navier-Stokes  limit  where  the  viscous  stress  and 
velocity  gradient  are  related  by  r  =  fidufdy,  one  may  define  a  new  spatial  variable 
ds  =  dy/jj,  and  write  r  =  dufds.  Because  t  is  a  constant  across  the  layer  in  Couette 
flow,  one  concludes  that  the  velocity  profile  must  be  linear  in  the  new  spatial  variable 
s.  However,  one  does  not  know  a  priori  whether  the  result  should  also  apply  to 
rarefied  flow.  The  figure  shows  that  straight  lines  are  in  fact  found  for  the  velocity 
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profile  for  all  degrees  of  rarefaction,  ranging  from  Maxwell  molecules  to  the  hard- 
sphere  limit.  In  addition,  it  shows  that  all  the  straight  lines  intersect  at  a  common 
point.  Because  the  Navier-Stokes  limit  is  represented  by  a  diagonal  line  connecting 
the  corners  of  the  square  box  enclosing  the  plot  and  the  free-molecule  limit  is  merely  a 
vertical  line,  and  its  position  can  be  predicted  theoretically,  the  intersection  point  can 
therefore  be  established  uniquely.  However,  a  most  remarkable  result  clearly  seen  in 
the  figure  is  the  fact  that  all  the  lines  pass  through  this  one  point.  The  temperature 
variable  is  also  shown  plotted  in  the  same  way  as  for  the  velocity,  and  again  all  the 
curves  pass  through  a  common  point,  but  a  slight  curvature  is  seen  in  the  temperature 
data  indicating  a  small  quadratic  term  is  present.  At  this  point  we  have  an  exciting 
result  but  we  do  not  have  a  suitable  theory  that  can  be  used  to  predict  it.  The  figure 
suggests  strongly  that  the  velocity  slip  and  temperature  slip  should  also  be  given  by 
a  simple  theory.  Some  work  has  been  carried  out  along  these  lines  by  Denery  but 
nothing  has  been  found  that  is  as  elegant  as  the  figure  suggests. 

Figure  16  shows  a  comparison  between  the  extended  Lees- Liu  solution  for  hard 
spheres  and  the  predictions  given  by  the  DSMC  method.  As  can  be  seen,  a  small 
difference  appears  in  the  transition  regime,  which  seems  to  increase  as  the  temperature 
ratio  between  the  two  plates  is  increased.  Otherwise,  the  two  give  the  same  predictions 
for  the  free-molecule  limit  and  for  the  Navier-Stokes  limit.  Relatively  large  values  of 
the  temperature  ratio  and  the  Mach  number  were  chosen  for  display  in  the  figure  in 
order  to  strain  the  comparison,  otherwise,  for  small  values  of  the  parameters  the  two 
agree  very  well.  The  cause  of  the  difference  for  large  values  of  the  parameters  is  most 
likely  associated  with  the  assumed  two-sided  Maxwellian  distribution  function  used 
in  the  Lees-Liu  method. 

Figure  17  shows  a  comparison  of  the  velocity  distribution  functions  given  by  the 
two  methods  and  it  is  clear  that  the  two-sided  character  of  the  assumed  theoretical 
distribution  is  not  fully  supported  by  the  DSMC  simulated  results.  The  general  trend 
is  certainly  present  but  the  details  are  significantly  different.  On  the  other  hand,  the 
Lees-Liu  velocity  distribution  function  was  originally  introduced  for  modest  degrees 
of  translational  nonequilibrium;  and  for  that  case,  the  discontinuity  seen  in  the  figure 
would  be  much  smaller  and  the  assumed  distribution  would  then  more  nearly  match 
the  DSMC  result.  A  complete  account  of  all  of  the  above  work  on  Couette  flow  can 
be  found  in  Denery’s  thesis  [4]. 
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PUBLICATIONS 


Three  student  theses  were  completed  during  the  period  of  the  grant.  These 
are  the  ones  listed  in  the  references  as  Fallavollita  [3],  Denery  [4]  and  Goswami  [6]. 
Likewise,  three  archive  publications  appeared  in  print  in  the  period,  and  these  are  the 
ones  by  Woronowicz  et  al.  [7],  BaganoIF  [8],  and  Fallavollita  et  al.  [9].  In  addition, 
an  AIAA  paper  by  Goswami  et  al.  [10]  also  appeared  in  the  same  period.  It  is  fully 
expected  that  two  or  more  journal  papers  will  soon  be  submitted  as  an  outgrowth  of 
Denery’s  work  and  likewise  several  should  also  result  from  the  work  of  both  Goswami 
and  Lou. 


STUDENTS  SUPPORTED 

A  total  of  three  graduate  students  were  supported,  at  various  levels,  in  the 
18-month  period  of  the  grant.  Terry  Denery  and  Tawei  Lou  received  their  principal 
support  from  the  grant  and  Craig  Dutweiller,  who  just  recently  joined  the  research 
group  to  pursue  a  Ph.D.  program,  received  partial  summer  support.  Douglas  Dahlby 
has  carried  the  principal  responsibility  for  continuing  the  development  of  our  parallel 
computer  code,  and  has  contributed  greatly  in  making  it  possible  to  carry  out  sim¬ 
ulations  having  many  practical  elements,  but  because  he  has  been  supported  by  an 
NSF  fellowship  he  does  not  appear  in  this  accounting. 
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Fig.  1.  Density  distribution  about  a  generic  wing-body  SSTO  vehicle  in  a  rarefied  flow, 
with  M  =  10,  Kn  =  0.02  and  a  =  20  degrees. 


Fig.  2.  Distribution  of  heating  flux  for  the  generic  SSTO  body  geometry  and  conditions 
considered  in  figure  1,  with  a  =  0. 


Fig.  3.  Velocity  field  as  seen  in  the  cross-section  plane  located  at  60%  of  the  wing  root 
chord.  Body  geometry  and  conditions  cure  the  same  as  those  for  figure  1,  except 
M  =  0.8. 


Fig.  4.  Velocity  field  as  seen  in  the  cross-section  plane  located  at  90%  of  the  wing  root 
chord.  Body  geometry  zuid  conditions  are  the  ssune  as  those  for  figure  1,  except 
M  =  0.8. 
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Fig.  5.  Effectiveness  of  restart  for  a  rarefied  Prandtl-Meyer  flow  having  a  Knudsen 
layer.  Time  development  of  a  conventional  start  (dashed  curve)  versus  a  restcirt 
(solid  curve)  is  shown  for  the  rms  temperature  difference  with  respect  to  a 
reference  solution. 


Fig.  6.  Effectiveness  of  restart  for  a  rarefied  supersonic  flow  past  a  circular  cylinder. 

Time  development  of  a  conventional  start  (dashed  curve)  versus  a  restart  (solid 
curve)  is  shown  for  the  rms  temperature  difference  with  respect  to  a  reference 
solution. 
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Fig.  7 


Location  of  the  largest  errors  in  the  early  stages  of  a  restart  simulation  for  the 
case  of  supersonic  flow  past  a  circular  cylinder.  The  quantity  displayed  is  the 
rms  density  error  between  the  restart  and  reference  solutions. 
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ect  of  a  low  number  of  particles  per  cell  in  a  simulation,  where  the  test  studies 
ished  curves)  made  use  of  1  particle  per  cell  in  the  free-stream  flow  with  (a) 
macrocells  and  (b)  macrocells,  while  the  reference  run  (solid  curves)  made 
!  of  16  particles  per  cell  and  no  macrocells.  The  variable  shown  is  temperature 
1  the  flow  is  about  a  circular  cylinder. 


Fig.  9.  Effect  of  spatial  resolution  (cell  size)  and  the  number  of  particles  used  in  a 
simulation  (particles  per  cell)  on  the  queility  of  the  simulated  results,  where 
improved  spatial  resolution  for  a  fixed  number  of  particles  is  shown  in  (a)  and 
increase  in  the  number  of  particles  holding  the  spatial  resolution  fixed  is  shown 
in  (b). 
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Fig.  11.  Comparison  of  a  hybrid  solution  (solid  curve)  with  a  full  DSMC  solution  (dashed 
curve)  at  the  time  the  reflected  shock  wave  just  passes  the  interface  position 
given  by  the  vertical  line. 


Y-Direction  Y-Direction 


Time  step  =  200  Time  step  =  500 
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Fig.  12.  Compaxison  of  a  hybrid  solution  (solid  curve)  with  a  full  DSMC  solution  (dashed 
curve)  for  the  problem  of  an  impulsively  started  flat  plate  amd  for  two  time 
periods:  am  eaurly  time  (200  time  steps)  when  r/p  <  0.1;  aind  a  later  time  (500 
time  steps)  when  r/p  >  0.1. 
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Normalized  Energy  Flux  Normalized  Shear  Stress 


Fig.  13.  Comparison  of  a  commonly  used  approximation  (dashed  curve)  with  the  full 
solution  (solid  curve)  for  Couette  flow,  M  =  0.8  and  the  hard-sphere  molecular 
model.  The  Navier-Stokes  result  (dotted  line)  is  included  for  reference. 
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Fig.  14.  Comparison  of  a  commonly  used  approximation  (dashed  curve)  with  the  full 
solution  (solid  curve)  for  Couette  flow,  M  =  2.0  and  the  hard-sphere  molecular 
model.  The  Navier-Stokes  result  (dotted  line)  is  included  for  reference. 
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Fig.  15.  Velocity  and  temperature  profiles,  as  determined  by  the  full  solution  for  Couette 
flow,  when  using  the  spatial  variable  defined  by  ds  =  dy/^.  The  upper  plate 
has  a  temperature  10  times  that  of  the  lower  plate  and  its  Mach  number  is  2.0, 
based  on  the  lower  plate  temperature. 


Fig.  16.  Comparison  of  the  full  solution  (solid  curve)  with  results  from  a  DSMC  simu¬ 
lation  (open  circles)  and  the  Navier-Stokes  prediction  (dashed  line).  The  case 
is  for  Couette  flow,  M  =  2.0,  and  the  hard-sphere  molecular  model. 
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Fig.  17.  Comparison  of  the  velocity  distribution  function  as  determined  by  a  DSMC 
simulation  with  that  found  from  the  full  solution  using  the  Lees-Liu  two-sided 
Meixwellian.  The  data  are  for  the  centred  plane  in  Couette  flow,  with  M  = 
2.0,  Kn  =  0.1  eind  a  temperature  ratio  of  10. 
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